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Abstract 

Magnetic resonance imaging (MRI) is a principal modality of modern medical imaging, which provides a wide 
spectrum of useful diagnostic contrasts, both anatomical and functional in nature. Like many alternative imaging 
modalities, however, some specific realizations of MRI offer a trade-off in terms of acquisition time, spatial/temporal 
resolution and signal-to-noise ratio (SNR). Thus, for instance, increasing the time efficiency of MRI often comes at the 
expense of reduced SNR. This, in turn, necessitates the use of post-processing tools for noise rejection, which makes 
image de-noising an indispensable component of computer assistance diagnosis. In the field of MRI, a multitude 
of image de-noising methods have been proposed hitherto. In this paper, the application of a particular class of 
de-noising algorithms - known as non-local mean (NLM) filters - is investigated. Such filters have been recently 
applied for MRI data enhancement and they have been shown to provide more accurate results as compared to many 
alternative de-noising algorithms. Unfortunately, virtually all existing methods for NLM filtering have been derived 
under the assumption of additive white Gaussian (AWG) noise contamination. Since this assumption is known to fail 
at low values of SNR, an alternative formulation of NLM filtering is required, which would take into consideration 
the correct Rician statistics of MRI noise. Accordingly, the contribution of the present paper is two-fold. First, it 
points out some principal disadvantages of the earlier methods of NLM filtering of MRI images and suggests means 
to rectify them. Second, the paper introduces a new similarity measure for NLM filtering of MRI Images, which is 
derived under bona fide statistical assumptions and results in more accurate reconstruction of MR scans as compared 
to alternative NLM approaches. Finally, the utility and viability of the proposed method is demonstrated through a 
series of numerical experiments using both in silico and in vivo MRI data. 
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I. Introduction 

Magnetic resonance imaging (MRI) is considered to be one of the most advanced modalities of modern medical 
imaging, which excels in providing a wide spectrum of useful diagnostic contrasts [1 1. Since the latter constitute an 
intensity-coded representation of biological properties of studied tissues/organs, the precision with which a contrast 
represents its associated biological property plays a decisive role in tissue characterization and early diagnosis. This 
fact establishes the value of post-processing techniques which aim at improving the signal-to-noise ratio (SNR) of 
diagnostic MR images, while preserving the integrity and consistency of their anatomical content. 

Unfortunately, in virtually all realizations of MRI, attaining higher spatial resolution entails using longer acqui- 
sition times. Apart from being highly undesirable from the perspective of patients' comfort and compliance, longer 
acquisition times lead to motion-related artifacts, which are the main foe of cardiac and diffusion MRI j2)-|5j. 
On the other hand, the reduction of acquisition time results in a loss of the spatial resolution as well as in an 
amplification of measurement noises. The latter tend to obscure and mask diagnostically relevant details of MR 
scans, thereby necessitating the application of efficient and reliable tools of image de-noising |6|. 

The current arsenal of image de-noising methods used in MRI is immense, which makes their fair classification 
a non-trivial task. For this reason, only three groups of de-noising methods which are germane to the present 
developments are mentioned below, while the reader is referred to the references therein for a more comprehensive 
literature review. In particular, the first group of de-noising algorithms for MRI encompasses variational methods, 
which are implemented through the solution of partial differential equations (PDE) |7)-fT0). Thus, for example, 
(7) suggests an adaptation of the classical anisotropic diffusion filter of [11] for noise reduction and enhancement 
of object boundaries in MRI. On the other hand, the de-noising method of JS) is based on minimization of an 
original cost functional, whose associated gradient flow has the form of a fourth-order PDE. In (9j, information 
from both the body coil image and surface coil image are incorporated in the form of data fidelity constraints. 
Finally, | [T0| introduces a maximum- a-posteriori (MAP) technique using a Rician noise model in combination with 
spatial regularization. 

A different group of de-noising methods takes advantage of the sparsifying properties of certain linear transforms 



p2|-|21 1 (for a comprehensive review of such methods, the reader is also referred to (22J). Thus, for instance, 
the method of fl6[ is based on wavelet thresholding applied to squared-amplitude MR images, supplemented by 
"unbiasing" of the scaling coefficients to account for the non-central chi-square distribution statistics. A different 



(robust) shrinkage scheme in the domain of a wavelet transform is proposed in |19|. Using a different line of 
arguments, the wavelet de-noising method of ]T8) is applied to complex-valued MR images. Finally, in pT) , the 
MR images are enhanced by means of a wavelet-domain bilateral filter. 
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A third group of image de-noising algorithms is based on the concept of non-local means (NLM) filtering, which 



was originally proposed in |23)-||25|, with its later improvements reported in |26|, p7| . As a general rule, NLM 
filters estimate a noise-free intensity of a given pixel (source pixel) as a weighted (linear) combination of the rest 
of the image pixels (target pixels). Here, the weights of the linear combination are determined based on a similarity 
measure (SM) between the neighbourhoods of the target and source pixels. As a result, the performance of an 
NLM filter is largely determined by the optimality of a chosen SM with respect to the properties of the image 
to be enhanced as well as those of the measurement noise. Thus, for example, under the conditions of additive 
white Gaussian (AWG) noise contamination, the above-referred NLM filters have been shown to outperform many 
variational and wavelet-based filters in terms of noise removal and the quality of edge preservation. 



Motivated by the success of NLM filtering in general image processing, the works in [28], |29| have extended 
the Gaussian-mode NLM filters to MR imagery. Additional reports on the subject also include J30)-J32"), where 
the filters are applied to MR images, followed by subtracting an estimation bias from the results thus obtained. 

Central to the main subject of the present work is the fact that the Rician statistics of MR images converges to 
Gaussian as the SNR of the images goes to infinity |3|. For sufficiently high values of SNR, therefore, applying the 
Gaussian-mode NLM filters seem to be well justified |3T) . However, for relatively low values of SNR (as in the 
case with, e.g., diffusion weighted imaging), the Gaussian model ceases to be legitimate, and as a result, the NLM 
weights optimal for the Gaussian setting become sub-optimal. This fact suggests a need for an SM which remains 
optimal for a wide range of the SNR values. Accordingly, deriving such an SM constitutes the main objective of the 



present paper. To this end, we start with an analysis of the NLM filters recently proposed in [33], |34| and underline 
some of their properties which should be avoided in the case of MR image de-noising. Subsequently, based on 
the results of f33") , p4| , we propose a new formulation of the SM and its associated weights, and demonstrate its 
usefulness and viability through a series of experiments using both in silico and in vivo MRI data. 

Table [IJ summarizes the main abbreviations and notations used in the paper, whose remainder is organized as 
follows. Section [II] provides some necessary details on the image formation model of MR images and their noise 
statistics. Sections [TIT] and [TV] describe a number of principal approaches to NLM filtering and point out some of 
their problematic aspects in relation to an MRI setting. A new SM and the closed-form expressions for its associated 
weights are derived in Section [V] while Section VI details a method for applying these proposed weights to the 



noisy pixels of data images. Section |VII| compares the performance of the proposed algorithm with that of some 
alternative methods using both in silico and in vivo MRI data. Finally, the main results and conclusions of the paper 



are recapitulated in Section VIII 



II. Image Formation Model and Noise Statistics 

MR images are standardly acquired in the Fourier domain (i.e., the fc-space), followed by the procedures of 
frequency demodulation and inverse transformation, which result in corresponding complex-valued images, whose 
magnitude is subsequently displayed Q. In this case, if the frequency-domain data is contaminated by zero-mean 
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TABLE I 

List of Notations and Abbreviations 



Notations and Abbreviations 


Meaning 


Formula (if applicable) 


NLM 


Non-local means 


- 


NCCS 


Non-central chi square 


- 


SM 


Similarity measure 


- 


SNL 


Similarity measure for NLM 




SSM 


Subtractive similarity measure 




RSM 


Rational similarity measure 




CSM 


Correlation similarity measure 




SNL v, fe 


Subtractive SM for NCCS distribution 






Rational SM for Rician distribution 






Subtractive CSM for NCCS distribution 






Rational CSM for Rician distribution 





AWG noise, the complex amplitude M of the noisy observation Aexp{ta} + N, with N — N r + iNi , is given by 

M = ^/(Acosa + N r ) 2 + (A sin a + iV,) 2 , (1) 

where A stands for the true image amplitude, while N r and JVj are mutually independent AWG noises of standard 
deviation a, and a £ [0, 2ir) is an arbitrary phase shift. In this case, M can be shown to follow the Rician conditional 
distribution model that is given bj[[] J3), Q 

J^exp{-4t ? !}/ (^), m>0 
PM\A(m\a) = S ( 2 ) 
0, otherwise. 



where J denotes the 0*' l -order modified Bessel function of the first kind. Figure 1(a) depicts several typical shapes 
of PM\A{ m \ a ) corresponding to a range of the values of A and a = 1. As can be seen from the figure, for A > 3a, 
the Rician probability density function starts closely resembling that of a Gaussian random variable [3|). However, 
for lower values of A, the density PM\A( m \ a ) becomes more asymmetric and protrudently heavy-tailed. Specifically, 
for A = 0, M follows a Rayleigh distribution model. 

The Rician nature of pm\a m Q renders impractical a straightforward application of many filtering strategies. 
This is because of the highly-nonlinear relation between the expectation £{M} of M and A. Specifically, 

£{M} = a^j2L 1/2 {-A 2 /2a 2 ), (3) 

where L v (x) denotes a Laguerre polynomial which, for v = 1/2, is given by 



Li,2(x) = e*/ 2 [(l-x)I Q (~)-xh (-|) 



(4) 



1 Here and hereafter, we use the standard statistical formalism for denoting random variables and their associated realizations by capital letters 
and their lower-case counterparts, respectively. 
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At the same time, a normalized version G = (M/a) 2 of the squared magnitude M 2 can be shown to be distributed 
according to a non-central chi square (NCCS) distribution with two degrees of freedom and parameter F = (A/a) 2 , 
whose conditional density is given by 

Pg\H9\J) = < (5) 
0, otherwise, 



where / 6 R + . Figure 1(b) shows a number of typical shapes of pQ\p corresponding to a set of different values 
of F. A better understanding of this figure can be derived from the fact that ([T| suggests 

G = F + 2VF£ + 77, (6) 

where £ := (iV r cosa + Ni sin a) /<r and T) := (N 2 + N 2 )/a 2 . Thus G can be viewed as a noisy version of F, 
where the noise has both additive and multiplicative components. Specifically, it should be noted that £ obeys a 
normal distribution with zero mean and unit variance, while 77 follows an exponential distribution with its mean 
and variance equal to 2 and 4, respectively. Moreover, the expectation of G now has a very simple relation to F, 
which is given by 

£{G} = F + 2. (7) 

It is the simplicity of |7]) which has been a principal impetus for the development of various de-noising methods, 
which have been applied to the squared magnitude G, rather than to its original value M. In one way or another, all 
these methods aim at recovering a close approximation of the average value £{G}, followed by the estimation of 
F through the subtraction of the global bias of 2. Note that, once an estimate of F has been obtained, its associated 
amplitude A in ([TJ can be recovered through taking the square root and re-normalization. These facts will be useful 
in the sections below. 



October 28, 2011 



DRAFT 



5 



III. Non-Local Means Filter 



The concept of NLM filtering was first proposed in [23| for the case of zero-mean AWG noise contamination. 
Let X s and Y s denote the intensities of the original (X) and observed (Y) image, respectively, corresponding to 



pixel s € I = {0, 1, . . . , N — 1}, where N denotes the total number of pixels in the image. Then, the filter of |23 1 
assumes Y s to be a realization of a mixture of ergodic processes, in which case the value of X s can be estimated 
by averaging {lt}tej,i where J s C I denotes the set of pixel indices whose associated intensities are distributed 
identically to Y s . It goes without saying that, in a real-life scenario, one might not be given an oracle who would 
provide us with all possible sets { J s } se j. In this case, it makes sense to replace the binary (hard) weighting by a 
fuzzy (soft) one, and compute an estimate X s of X s according to 

X s = W M Y t, with C s = ^ w„, t , (8) 

s tei tei 

where w s ,t > quantifies the "contribution" of a target pixel t E I to the estimate of the source pixel s. Note that, 
ideally, the weights w Si t should reflect a degree of similarity between the image intensities in the vicinity of pixels 
s and t of the original image. Thus, for example, the original choice of (23) was 



w Stt = exp 



(-tE&I^-*-^-*! 2 ), ( 9 ) 

V ken / 



where ft is the index set representing a symmetric neighbourhood of the centre of image coordinates. (Thus, for 
example, with s corresponding to two spatial coordinates s ■<-» (x,y), il could be defined as ft = {(x,y) : \x\ < 
L x , \y\ < L y } for some positive integers L x and L y .) The "fine-tuning" parameters {/3fc}fc e o in Q are intended 
to weight the domain of summation and they should be chosen to satisfy J^ken @ k ~ ^' wn il e h > controls 
the overall amount of smoothing imposed by the filter. Specifically, the higher values of h tend to result in overly 
smoothed output images, whereas smaller values produce rather milder filtering effects. As a general rule, an optimal 
value of h should be chosen adaptively according to the level of noise in Y. 

To facilitate our considerations, we note that the NLM weights in |9]l can be alternately expressed as [33) , [34) 

w,, t = n( SNL -A*)^. 
with SNL S t k being a Gaussian SM defined as 

SNL S)t , fe = exp (-|n_fc - Y t _ k \ 2 ) . (11) 

It is important to note that the value of SNLg^fe is always bounded between and 1, and it can be shown (see 
||33j, (35) for more details) that its choice in (jTOJ and ( fTT) is optimal in the case of AWG noise contamination. 

It should be noted that the use of weights {w s t } for linearly weighting the noise samples Y t is not the only way 
in which the weights can be "intermingled" with noisy data to yield an NLM estimate X s similar to ([8j. Thus, for 
instance, an alternative way of weighting is used in [33], [34- 1 based on the maximum likelihood (ML) framework. 



Apart from the ML-based weighting scheme, the works in (33), (34 1 suggest a unified approach to computation of 
the optimal weights based on the formal (and, in general, non-Gaussian) statistical properties of the original image 
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as well as of measurement noise. Despite the generality of the above formulation, however, its application to the 
case of MR imagery results in SMs which possess a number of undesirable properties. In the section that follows, 
these properties are brought under consideration, followed by the derivation of an original methodology that allows 
one to fix them. 



IV. Statistical approaches to computation of SM for MRI 



In 1 33 1, it was suggested to set the SM SNLg^^ to be equal to the posterior probability of X s - k = X t -k 
conditioned on observations of Y s - k and Y t -k- Formally, 

SNL s , t , fc = P(X s _ k = X t _ k \Y s _ k ,Y t _ k ), (12) 

which, for the case (3 k /h — 1, Vfc G ft, leads to following definition of the NLM weights 

w>,t = II P{X s -k = X t _ k \Y s _ k ,Y t _ k ). (13) 
ken 

It is worthwhile noting that, under the assumption of statistical independence of the intensities of the original image 
X, the weights in {[3) can be viewed as the posterior probability of the image patches {^ s -fc}feeo and {^t-k}ken 
to consist of the same intensities, conditioned on the observation of their corresponding noisy values {Y s ^ k }ken and 
{Y t ~k}ken, respectively [33]. Although the assumption of statistical independence is an obvious oversimplification, 
it is often employed in NLM filtering to render the final estimation scheme computationally feasible. 

It should be noted that the SM in ( fl2] > seems to have a serious theoretical flaw in the case of continuous 



random variables X s _ k and X t - k , in which case the SM is always equal to zero |36 p. 111]. To overcome this 
difficulty, it was suggested in ]34) to introduce an auxiliary random variable U k = X s - k — X t -k and set the SM 
SNL S t k to the value of the conditional density Pu k \Y 3 - k ,Y t - k i u k\Us-k, Ut-k) at u k — 0. Alternatively, one can 
use a different auxiliary variable V k = X s - k j X t - k , in which case the SM can be set to be equal to the value of 
Pv k \Y s k .Y t k (vk\ys-k,yt-k) at v k = 1 0341 . For the convenience of referencing, the above SMs will be referred 
below to as the subtractive and the rational SMs, respectively. Note that, although similar in their underlying 
philosophy, these measures lead to substantially different de-noising schemes, as it is detailed below. 

In the present work, we explore both the subtractive and rational SMs for two different types of input data, 
namely for the measured magnitude M and its squared normalized version G. The main contribution of the next 
propositions is to provide closed-form expressions for the SMs which result from using a subtractive U k on G 
images, and a rational V k on M images. Unfortunately, for the remaining two combinations (viz., "rational" G 
and "subtractive" M) it does not seem to be possible to derive closed-form expressions as well. In these cases, 
the measures need to be computed numerically - the approach which should be avoided in practice due to its low 
computational efficiency. 

Proposition 4.1: Let G = (M/a) 2 be a squared and normalized version of the MR magnitude image M as given 
by ([TJ. Moreover, let F = (A/a) 2 , where A denotes the true signal amplitude. Then, the subtractive SM SNL^ 1 ] k 
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is given by 



SNL« k = p Fs _ k - Ft _ h \ G3 _ k , Gt _ k (0 | g s - k ,gt-k) = \ e^-- ' 



(14) 



The proof of Proposition 4. 1 is provided in Appendix |Aj while Fig. 
corresponding to various values of gt-k and some fixed values of g s - k - 



2(a) 



(i) 



shows a number of SNL;, '. , curves 



Tfl 




20 25 30 



(a) 



(b) 



Fig. 2. (a) Subtractive similarity measure SNL^ , ; (b) Rational similarity measure SNL' 2 ? fc . 



This effect is due to the 



Observing Fig. |2(a)| a number of critical remarks are in order. 

1) The smaller the value of g s -k is, the narrower is the effective support of SNL^t. 
signal-dependent Gaussian noise 2\/F£ in (|6j. Indeed, smaller values of g s -k suggest smaller values of their 
related f s -k, and therefore smaller values of the above mentioned Gaussian noise component. In such a case, 
SNL^ k becomes more sensitive to the value of g s -k ~ gt-k (fast decay), which indicates that the difference 
fs-k ~ ft-k is likely to be small as well. On the other hand, the contribution of the (signal-dependent) 
Gaussian noise component becomes stronger for relatively large values of g s - k (and hence of f s -k)- m this 
case, SNLg£ fe has a slower convergence rate, since larger values of g s -k ~ gt-k no longer imply a larger 
discrepancy between f s -u and ft- k - 

2) The smaller the value of g s - k is, the more heavy-tailed is the behaviour of SNL^ fc . This effect can be 
attributed to the dominance of the exponential noise component r/ in (|6j. Moreover, when g s - k is large 
(which implies the dominance of the Gaussian noise), the SNL^ k curves appear to be more symmetric. 

The proposition that follows extends the previous results to the case of magnitude signals and a rational SM. 
Proposition 4.2: Let M be the MR magnitude image as given by ([TJ. Moreover, let A be the original signal 

(2) 



amplitude. Then, the rational SM SNL^ ; k is given by 



SNLg ifc =p A . 



. k /A t _ k \M s _ k ,M t _ k 



(1 | m s - k) m t - k ) = ^-1^ e-rt-'+t-Vlo (^^) • (15) 
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The proof of Proposition 4.2 can be found in Appendix |b| The plot of SNlA 2 : k is shown in Fig. 



2(b) 



where each 



curve is drawn with a fixed m s -k and varying m t -k, while a is set to be equal to 1. 
The main observations in the above case are: 

1) As the value of m s -k increases, the noise distribution tends towards Gaussian, and as a result, the similarity 

(2) 

measure SNL;, ; . becomes more symmetric and appears to have the shape of a Gaussian SM. 

(2) 

2) For relatively small values of m s _fc, on the other hand, the shape of SNL;, ' t fc is noticeably asymmetric, being 
heavy-tailed towards the right. 

Although SNLg j k and SNL^ k appear to reflect the main properties of their corresponding noise distributions, 
they share a number of critical drawbacks which we point out below. 

1) Neither SNL^ fe nor SNL^ k (as given by ( fB) and (15) , respectively) attain their maxima values at the 
point when the two arguments received by the measures are equal. Particularly, for a fixed value of g s -k 
(resp. m s -k), the SNL^ fc measure (resp. the SNL^ fe measure) is maximal at some g t -u < g s -k (resp. 
m t -k < m 3 -k). 

2) The maximal values SNL^ k and SNL^ 2 ] fc can attain depend on the values of the "source" intensities g s -k 
and m s _fc, respectively. In other words, the value of SNL^ 1 ] k (a, a) (resp. SNL* 2 ^ fc (/3, /?)) depends on the 
value of a (resp. j3). From a purely applicational point of view, this fact suggests that the measures are not 
scale invariant, and as a result, the weights w s ,t in are defined not only by how dissimilar compared 
intensities are, but also by their absolute values. 

(2) 

SNL^ ' t k can become unbounded which is not a favourable property of the 



2(b) 



3) As can be seen from Figure 
similarity measure. 

/i ^ /'n\ 

where 



Some of the above mentioned limitations of SNL^ 1 ] k and SNL^ 2 ] k have already been pointed out in 1 34 



the authors apply the method of J33) to multiplicative noise. The main conclusion which one can immediately draw 
from the discussion above as well as based on 134) is that neither SNL^ . nor SNL^ , is optimal to deal with 
the cases of Rician and/or NCCS noises, which are the most relevant types of noises in MRI. In the next section, 
we propose a new SM, which is free of the limitations of SNL^ k and SNL^ k mentioned above. 

V. Proposed Approach 

A. Subtractive SM for NCCS noise 

To derive the proposed SM in a consistent and intuitive way, we start with the case of the subtractive SM (SSM), 
which is defined as |34| 

SSM = PF B _ k -F t - k \G s - k ,G t - k (° I 9s-k,9t-k) = 

PF 3 - h \G 3 - k ,G t - k {f I 9s-k,gt-k)PF t - k \G s - k ,G t - k (f\9s-k,gt-k)df. (16) 
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For the case of a NCCS distribution, one alternatively has (see Appendix [A] for the details of derivation) 

SSM = SNL« =p Fs _ k -F t _ k]Ga - k ,G t - k (V I 9 s -k,9t-k) = 



PG s - k \F s - k (9s-k I f)PG t - k \F t _ k (9t-k I f)df. (17) 

In the above equation, it has been assumed that the prior probability PF z _ h (where z £ {s, t] and k g fl) is uniform. 
Moreover, it has also been assumed that, given the noisy intensity value at a particular location, its original (i.e., 
noise-free) value is conditionally independent of noisy intensities at different locations (viz., Pf, _ h \G xl -t,,G, 3 -k = 
PF. 1 _ fc |G 21 „ fc )- Finally, it also deserves noting that, for a general noise distribution, it can be shown that the relation 
in ( fPTj ) holds with the proportionality rather than the equality sign J33) , J34|. Specifically, 



SSM = p Fa _ k _ Ft _ klGs _ k , Gt _ k {0 | g s ^ k ,g t „ k ) oc / p G3 _ k \ Fa _ k (g s -k \ f)PG t _ k \F t _ k (9t-k I /) df. 

Jo 



(18) 



To overcome the limitations of SNL^ , and SNL^ > as per the discussion in the previous section, we interpret 
the right-hand side of ( fT8] l as an inner product between the likelihood functions L gsk (-) and L gtk (-), with L g (f) 
given by 

L g {f)=p G]F {g\f). (19) 

In other words, the SM derived from the probabilistic point of view has the form of the inner product of likelihood 
functions, where each of the likelihood functions is indexed by its corresponding noisy observation. Note that, in 
general, the likelihood functions L 9a k {-) and L gt k (-) have unequal norms, and as a result their inner product 
is not maximized when g s -k = 9t-k (which would be a natural and desirable property of an SM to have). To 
overcome this shortcoming, we suggest to normalize the inner product, thereby converting it into a correlation 
similarity measure (CSM) according to 



(, ^= ,,r £ Vifr^ii ' (2o) 



where (x, y) = J °° x(f)y(f) df and ||x||2 = \j (x, x). It is interesting to observe that CSM^ k can be also viewed 
as an inner product of two functions lying on the unit sphere in L 2 (1R + ). Therefore, CSM S t ^ has an interpretation 
of the cosine of the angle between the two functions. 

The CSM in (|20b turns out to be particularly advantageous in the case of MR imagery. First, it is free of all 



the major limitations of SNL^ 1 ] k and SNL^ 2 t ' fe as previously discussed. In particular, CSM s t .fc is always smaller 
or equal to 1, and it achieves its maximum value when g s ~k = 9t-h- Secondly, for the case of NCCS noises, the 
CSM measure can be shown to have a neat closed-form expression which is given by 



c\tt( 3 ) A) {y/9s-kgt-k) , 

SNL »,t.fc = ~/ r /- w where 3 = . g/2. (21) 

V h\9s-k) h\9t-k 



(3) 

A number of graphs of SNL^, f fc are shown in Fig. 3(a) where each curve is drawn with a fixed g s -k and varying 



9t-h- It can be seen from the graphs that the shape of each curve is similar to those in Fig. 2(a) However, unlike 



the plots in Fig. 2(a) each curve in Fig. 3(a) is maximized when gt-u = g s -k and the maximum value is the same 
for all the curves and is equal to 1. 
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Qt-k 



20 25 30 
Wf-6 



(a) 



(b) 



Fig. 3. (a) Proposed CSM SNL^ fc for the case of NCCS noise distribution; (b) Proposed CSM SNL^ fc for the case of Rician noise 
distribution. 



B. Rational SM for Rician noise 

To derive an expression for a CSM for the case of Rician noise, we first recall that the rational SM (RSM) is 
given by (34) 



RSM = p As _ k/At _ klMs _ kyMt _ k (l | m s _ k ,m t ^ k ) 

aPA s _ k \M s _ k ,M t - k ( a I m s-k,rn t _ k )p At _ k \ Ms _ k . Mt _ k (a \ m s ^ kl m t -k) da. (22) 
Using a method similar to the one discussed in the previous subsection, one can obtain 

/>oo 

RSM= / ap MB _ k \A a _ h {m s -k \ a) p Mt _ k \ At _ k (m t - k \ a)da, (23) 
Jo 

where the equality sign can be replaced by a proportionality sign for general noise distributions, as it was done 
in the previous subsection. Moreover, similar to the case with SNL^ fc , the integral in ( |23~| can be interpreted 
as a weighted inner product (x,y) a = x{a)y{a) a da , where a da can be viewed as a "modified" integration 
measure. Using this notation, we have 

RSM=(L ms _ h ,L mt _ k ) a . (24) 



As the final step, one can normalize the inner product in ( |24[ > in a way similar to (|20j), which leads to a different 
CSM for Rician noises, which is given by 

T C 

(4) _ J {- 



SNL 



2<r 2 



s,t,k 



(25) 



A number of SNL^ t k curves are shown in Fig. 3(b) for different (fixed) values of m s - k and a range of varying 



mt-k- Once again, one can observe that the curves are similar in shape as those of Fig. 2(b) However, unlike Fig. 
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2(b) 



each SNL^ 4 ] k curve is maximized when mt-k = m s-k an d the maximum value is equal to 1 in each case. 
There are two important facts about SNL^ k and SNL^ fc that deserve to be paid special attention. In particular: 
1) The values of SNL^ fe and SNL^ k (as given by ( JiTj ) and (25 i, respectively) are equal under the substitution 

\2 T« t~tt Kar M7/-vi-/-lc fha i/ilno rvF QT\TT 



9z-k — (mj-d/ff) . In other words, the value of SNL;, t fc for some arbitrary g s t k and (j 1 ^ will be equal to the 
value of SNL^ fc for m s ^ and m t k, whenever g s ^ = (™ s ,fe/cr) 2 and g t j, — (m f fc/cr) 2 , which is precisely 
the case at hand. In other words, the proposed SMs are equivalent under reparametrization G = (M/a) 2 . 
This fact suggests that, in terms of SNL^ fe and SNL^ k values, two patches of an original MR image and 
their corresponding squared and normalized versions are equally similar. 
2) Using the fact that, for sufficiently large x, it holds that 

I (x) w (26) 
V Zttx 



and plugging d26b into ( 25 i instead of the original Bessel functions, we obtain 



SNL « k « exp ( - K -V_f f - fc|2 | . ( 27 ) 



4(7 2 

The above approximation holds with a high precision for relatively large values of SNR (i.e., for ^ 3> 1). 
This is an exceptional property of SNL^ 4 t . , since it is known that Rician noise in MRI converges to Gaussian 
noise, when SNR increases. In such a case, the proposed SNL^ k measure converges to the form of ( fTT) , 
whose optimality for Gaussian noises was proven in p3|. 



VI. Bias Removal 

The original NLM filter of [23 1 employs the adaptive local averaging scheme of ([8]), in which the output intensity is 
computed as a weighted linear combination of the measured intensities, while the weights are determined adaptively, 
based on an SM in use. Unfortunately, the filter in (|8]l has some optimality properties in the case of additive Gaussian 
noises alone, while the efficiency of the filter is known to deteriorate for some more general types of measurement 
noise. In the latter case, one can adopt the more general filtering procedure that was used in (33], [34) , according 
to which an estimated value X s of some unknown X s is given by 

X s = axgmax^ w Sjt log p(Y t | X s ), (28) 
3 tez 



where Y t denotes the measured intensities and {w s j} are some predefined weights. The estimate in ( |28] l has been 
motivated by the work in (37) where it has been referred to as a weighted maximum likelihood (WML) estimate. 

Unfortunately, in the case when Y t in (28} follows either a Rician or an NCCS distribution, (28} does not seem 
to have a close-form analytical form, which necessitates the use of numerical optimization. Alternatively, one can 
take advantage of the relation in (7} to estimate £{G S } as 

£{Gs}~Yl Ws ^9u (29) 
tei 
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followed by estimating the related original amplitude A s = <jy/£{G s ~} — 2 as 

n 1/2 



A., 



iW. (30) 



where the max(-) operator is used to avoid complex estimates. Note that, even though applying the max(-) operator 
may seem like a purely ad-hoc procedure, it has been shown to be optimal in the ML sense in [38]. This estimate 
has also been used previously in pT) in the context of non-local means. 

Yet another option to estimate A s is to apply averaging to M (as opposed to G). Although the square of the 
average value is, in general, not equal to the average of the squared values, a common method of estimating A s is 
1301 

1/2 

= Af\ (31) 



-4., 



2 ^ 

-2a 2 , 



w s ,t m t 

tei 

Despite the fact that Etej w s.t m t] 2 l& 2 does not seem to be a legitimate estimate of £{G S }, the estimate A^ 
in ( (3T| often produce more accurate reconstruction results in terms of the mean-squared error (MSE) as compared 
to the estimate A^ in d30l>. To understand the reasons which underpin the above phenomenon, it is instructive 



to consider the following numerical experiment. Fig. 4(a) and Fig. 4(b) show the histograms of the estimates in 



( |3"0"1 > and ( pi") , respectively, which have been computed without applying the square root and the max(-) operator. 
With a slight abuse of notations, these estimates are referred to below as (Ai 1 ^) 2 and (A^) 2 (which, despite the 
square in their superscripts, are allowed to have an arbitrary sign). In both cases, the original amplitude A was set 
to be equal to zero (i.e. A = 0), the noise variance a was normalized to have the value of a = 1, the weights w s _ t 
were chosen to correspond to uniformly averaging 25 random realizations of g± and mt, respectiveljj^] and both 
histograms have been computed based on the results of 10 6 independent trials. Observing Fig. [4] one can see that 
both estimates provide outputs in the vicinity of the true squared amplitude A 2 = 0. However, while the histogram 
of (A^) 2 is shaped more or less symmetrically around the origin, the histogram of (A^) 2 is noticeably biased to 
the left, which suggests that the estimate in pTj ) tends to underestimate the true squared magnitude of A 2 = 0. In 
terms of the probabilities P^Ai 1 ^) 2 < 0) and P((A^) 2 < 0) with which the estimates yield non-positive outputs, 
one can therefore conclude that P{(A^) 2 < 0) < P{{A^) 2 < 0). (Thus, for example, in the case of Fig. Q 
P((ii 1) ) 2 < 0) fa 0.54, while P{{Af ] ) 2 < 0) 0.89.) This fact, in turn, implies that, after applying the max(-) 
operator, the estimation in pTj ) is much more likely to produce the exact reconstruction of A — as compared to 
the estimation in ( f30] >. 

In MRI, zero-valued amplitudes A are predominant at the background areas of MR images, which are normally 
devoid of water content. At such areas, therefore, the estimation in pO) should be expected to produce larger values 
of MSE as compared to the case of pTj ). We will have more to add to the subject in the following sections of the 
paper. 

2 In practice, 25 appears to be a typical size of the neighbourhood / used for NLM filtering. 



October 28, 2011 



DRAFT 



13 




VII. Results 

A. Reference Methods 

The performance of the proposed method has been compared with several standard and established algorithms. 



Specifically, as the first reference method, the total-variation filter of 1 39 1, implemented by means of the fast fixed- 
point algorithm of ]40) , has been used. The algorithm (which is referred to below as TVDN, an acronym of total 
variation denoising) has been applied to the magnitude MR data, followed by the bias-removal procedure specified 
by ( (3T| . As the second reference method, the wavelet-based method of [19| has been usecQ In what follows, this 
method is referred to as wavelet de-noising (WDN). As the final method used for numerical comparison, the NLM 
filter of p0| has been employed. Since the filter uses Gaussian weights to compute the similarity measure, we refer 
to this method as GNLM, an acronym for Gaussian NLM. In the case of all reference methods under comparison, 
their respective parameters have been set based on the guidelines specified in their associated papers. 

Two new approaches to NLM filtering of MR images are proposed in this paper. Specifically, the first approach is 
applied to squared-magnitude MR data using the subtractive CSM of ( pi) , followed by the bias-correction procedure 
given in ( f30] >. For the convenience of referencing, this filtering approach is referred below to as NLMS (with "S" 
standing for "subtractive"). The second approach, on the other hand, is applied to magnitude data using the rational 
CSM of ( p5| , followed by the bias -correction procedure of ( |3T| . In what follows, this filtering approach is referred 
to as NLMR (with "R" standing for "rational"). All the acronyms of the proposed and reference algorithms are 
summarized in Table |ll] 



3 The method was implemented using the code available at the author's webpage at 

Sof tware/MRIprogram . zip . 



http://telin. ugent . be/ " san ja/San ja_f iles / 
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TABLE II 

Acronyms of the proposed and reference algorithms 



Algorithm name 


Reference 


Input image type 


TVDN 




39 




M 


WDN 




19 




G 


GNLM 




30 




M 


NLMS 


Proposed 


G 


NLMR 


Proposed 


M 



B. Simulated Data 

In this subsection, the MRI images from the built-in MRI dataset available in the MATLAB® toolbox have been 
used as test subjects. Specifically, the denosing algorithms have been tested using the axial slices number 4, 7 and 



16 (shown in Figures 5(a) 5(b) and 5(c) respectively), which represent a spectrum of different cerebral structures. 
For quantitative comparison, simulated data have been obtained by subjecting the original test images to various 
levels of Rician noise. 




(a) (b) (c) 

Fig. 5. (a) Test slice #4, (b) Test slice #7 and (c) Test slice #16 of the MATLAB® MRI database. 



1 ) Data Generation: For the current test setup, given a noise-free intensity A k of the original image A, its noisy 
counterpart Mfe can be simulated according to 

M k = ^{A k +n r y + n 2 t (32) 

where n r , m ~ Af(0,a 2 ) are independent Gaussian random variables, which are also assumed to be independent 
across the image domain. In our simulation study, the standard deviation a was set to 10 percent of the maximum 
value of the original images, thereby resulting in a (peak) signal-to-noise ratio of SNR = A max /a = 10. 

2) Performance Metrics: Evaluating the quality of denoising in medical imaging is of subjective nature, as it is 
often based on the particular requirements of a medical expert. For this reason, in many cases a noisy version of an 
image is given preference over a de-noised version, as denoising procedures have an inherent risk of removing small 
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structures from images. Hence, it is difficult to define an objective criterion for evaluating the quality of medical 
images. Nevertheless, there is a number of standard evaluation metrics used in the literature, some of which we 
adopt in the present study. Specifically, one of such metrics is the root mean square error (RMSE), which can be 
expressed in dB as given by 



RMSE = 20 log! 



N 

2 



k=i 



1/2 



(33) 



where = A^ — A^ denotes the difference between the original intensity A^ and its estimated value A^ at position 
k, and N stands for the total number of image pixels. 

When estimate Ak is biased, the mean value of may not be equal to zero, in general. In this case, it makes 
sense to replace in ( [33) ) by its centred version — e, with e being the sample mean of given by 



N 

1 

e 



1 N 



N 
k=i 

The resulting metric is called the centred RMSE (cRMSE), and it is formally defined as 

1/2 



cRMSE = 20 log 



10 



1 E 



N 

2 



N ^ ^ ~ 6 I 

k=l 



(35) 



It should be noted that, while structurally similar, the RMSE and cRMSE metrics provide different quantitative 
assessment in the case of biased estimation. Consequently, the analysis and comparison of both these metrics can 



be helpful in evaluating the performance of the de-biasing procedures detailed in Section VI 

Despite the relatively straightforward interpretation offered by the RMSE and cRMSE metrics in ( [33) ) and ( f3"5] l, it 
has been argued in pTj that these metrics may not always adequately represent the effect/size of certain estimation 
artifacts/noises. As an alternative, a different performance metric - called the structural similarity index (SSIM) - 
has been proposed. This metric has also been used in our comparative study. 

3) Details on the choice of parameters: In order to facilitate the reproducibility of the proposed algorithms and 
their results, some principal details on the choice of algorithms' parameters are specified next. In particular, all the 
NLM algorithms have been implemented using a square 5x5 neighbourhood £1 and an 11 x 11 search window. 
(Note that the latter suggests that the averaging in |8]l was carried out only over pixels t which were located within 
an 11 x 11 window centred at pixel s.) The weights j3k in ( fTO) were defined to correspond to a separable binomial 
mask of the third ordei|^] while the tuning parameter h was set to be equal to 0.4. Note that, in general, it has been 
found empirically that the proposed algorithms performed the best for h 6 [1/3, 1/2]. 



4) Comparative analysis of algorithm performance: Table III summarizes the values of the performance metrics 



of Section VII-B2 obtained for various tested images and different denoising algorithms (with the best results being 
accentuated with bold characters). Furthermore, for the sake of visual comparison, the image restoration results are 
shown in Fig. [6] Fig. [7] and Fig. [8] which have the same composition. Namely, Subplots (a) of the figures show 

In practice, this mask can be computed as vv , with v = 2" 4 ■ [1 4 6 4 1] T . 
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TABLE III 

Performance metrics of the result of denoising using different algorithms 



Image 


MRI Slice 4 


MRI Slice 7 


MRI Slice 16 




RMSE 


cRMSE 


SSIM 


RMSE 


cRMSE 


SSIM 


RMSE 


cRMSE 


SSIM 


Noisy 


21.06 


18.6175 


0.4187 


20.8266 


18.6651 


0.4193 


20.8274 


18.7672 


0.3732 


TVDN 


13.8097 


13. 8094 


0.9028 


13.1824 


13.1777 


0.8974 


13.4104 


13.4078 


0.8878 


WDN 


16.0125 


14.9681 


0.4849 


15.3424 


14.4047 


0.5158 


14.9729 


13.9367 


0.4890 


GNLM 


13.4437 


13.4261 


0. 9128 


12.4910 


12.4715 


0.9104 


11.3428 


11.3431 


0.8920 


NLMS 


13.3571 


13.1606 


0.7066 


12.5003 


12.3805 


0.7543 


11.9151 


11.6872 


0.7163 


NLMR 


12.6287 


12.6115 


0.9270 


11.8773 


11.8554 


0.9204 


10.6613 


10.6534 


0.9110 



the noisy data images, while Subplots (b)-(f) show the recovered images obtained using the TVDN, WDN, GNLM, 
NLMS and NLMR algorithms, respectively. 



As can be seen from Table III as well as from its related figures, the NLM algorithms produce better reconstruction 
results as compared to TVDN and WDN. Moreover, among the NLM algorithms, NLMR provides better performance 
than NLMS and GNLM in terms of all the performance measures. It is worthwhile noting that the performance of 
NLMS is comparable to that of NLMR and GNLM in terms of the RMSE and cRMSE measures, while the former 
presents significantly lower values of SSIM as compared to the other two. To further explore this phenomenon, 
Fig. [9] shows the SSIM maps pT) computed for the images of Fig. [6] The SSIM maps represent the local values 
of SSIM (with their brighter intensities indicating stronger resemblance between the reconstructed and original 
MRI images), and hence they are particularly suitable for analyzing the spatial distribution of reconstruction errors. 
Thus, for example, Subplot (a) of Fig. [9] corresponds to the noisy data image, in which case the SSIM values 
at the image background are negligibly small, indicating little resemblance between the original (monotone) and 
measured (noisy) data. Similarly, the SSIM map corresponding to the NLMS reconstruction (as shown in Subplot 
(e) of Fig. [6]i has considerably darker background values in comparison to the GNLM and NLMR maps (shown 
in Subplots (d) and (f), respectively). The NLMS reconstruction's reduced background resemblance is mainly due 
to the bias subtraction (thresholding) procedure discussed in Section |VI] Even through a visual inspection fails to 
find a difference between the original and the reconstructed background, the accumulation of small errors over a 
relatively large number of background pixels adversely affects the average value of the SSIM metric, as indicated 



by Table III 



An additional important observation can be made through comparing the values of RMSE and cRMSE metrics 
obtained using different reconstruction algorithms. In particular, we first note that the values of RMSE and cRMSE 
corresponding to the noisy data are not identical - the fact which indicates the presence of a non-zero bias in 
the measurement noise. At the same time, these metrics have approximately equal values for the case of TVDN, 
GNLM, NLMS and NLMR reconstructions, while being noticeably different in the case of WDN. This fact suggests 
that the latter method is inefficient in removing the constant bias in the reconstruction error. 
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(d) (e) (f) 

Fig. 6. (a) Noisy version of slice #4; (b)-(f) Reconstruction results obtained using TVDN, WDN, GNLM, NLMS and NLMR, respectively. 



C. Experiments with real-life data 

Reconstruction of real-life MRI images has been the next step in our comparative study. To this end, the data set 
of |19 | have been used herein. The data were obtained at the University Hospital of Ghent and it is publicly available 
at |http : //telin . ugent . be / ~san ja/ San ja_f iles/ Sof tware/MRIprogram. zip| The data contains 
a sagittal and an axial scan of a human brain, which are shown in Fig. |10(a)| and Fig. |1 l(a)| respectively. 

The reconstruction results obtained for each of the tested images using the proposed and reference methods are 



shown in Subplots (b)-(f) of Fig. 10 and Fig. 11 respectively. From these figures, it can be seen that the proposed 
algorithms result in higher-contrast reconstructions of better visual clarity as compared to the reference approaches. 
The difference is particularly evident for the case of Fig. [TT] where the proposed algorithms result in less noisy 
images, while exhibiting higher effective resolution and contrast. 

VIII. Conclusion 

The present paper has proposed two novel NLM-based methods for enhancement of MR images. More specifically, 
the paper introduced a new definition of the NLM weights, which takes into consideration the true Rician statistics 
of measurement noises. The proposed definition has been shown to overcome some unfavourable characteristics of 
the weights previously proposed in the MRI denoising literature. Subsequently, the current paper provided closed 
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(a) 



(b) 



(c) 




(d) 



(e) 



(f) 



Fig. 7. (a) Noisy version of slice #7; (b)-(f) Reconstruction results obtained using TVDN, WDN, GNLM, NLMS and NLMR, respectively. 

form expressions for the weights corresponding to both subtractive and rational similarity measures. The utility of 
the proposed algorithm has been demonstrated through a series of computer simulations and real-life experiments. 
Based on the obtained results, one can find that the proposed algorithm has been able to provide better reconstruction 
results as compared to a number of established reference approaches. 



Appendix A 

Subtractive SM for non-central chi square statistics 
We are interested to evaluate SSM = PF s _ k -F t _ k \G a - k ,G t - k (0 I 9s-k> 9t-k)> which is defined as |34| 

p oo 

SSM= / PF,- h \G a - h ,G t -M\9s-k,9t-k)VF t - h \G.- h ,G t -M\9s-k,9t-k)df, (36) 
Jo 

where it has been assumed that F s _fc and Ft-k are conditionally independent, given the noisy data G s -k and Gt-k- 
Moreover, assuming 



PF» 1 - k \G :s:L - k ,G* 2 - k (f\9z 1 -k,gz 2 -k) =PF zl - k \G n - h {f\9z 1 -k), z l7 z 2 G {s,t} 
and using Bayes theorem, one obtains 

Jq 00 PF.-k(f) PF t - k (f) PG 3 - k \F a - k (g.s-k | f)PG t - k \F t - k (gt-k I f)df 



SSM = 



PG s - k (g s -k) PG t - k (gt-k) 



(37) 



(38) 
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(a) 



(d) 



(b) 



(e) 



(c) 




(f) 



Fig. 8. (a) Noisy version of slice #16; (b)-(f) Reconstruction results obtained using TVDN, WDN, GNLM, NLMS and NLMR, respectively. 



If no prior information about the original intensities is known, one can assume that F s -k and Ft-k follow a uniform 
distribution |33| , |34) . Also, for the case of an NCCS distribution, PG*_ fe (<?z-fc)> where z € {s,t}, can be written 

as 



PGz- k {9z-k) = / PG*- k \F z - k (9z-k I fz-k)PF z _ k (f)df = 
JO 

l-oo i 



(39) 



= PF z - k (f), 

where the second line in ([39) follows from the uniform density assumption of F z -k and the integration result is 
obtained by noting that PG x _ k \F z - h can a l so be looked upon as a probability density function of F z -k with g z -k 
as its parameter. Hence, in the case of NCCS noise distribution, the expression for SSM is given by 

(i) . 



SSM = SNL 



Substituting y = \ff results in 



,t,k 



PG s - k \F s - k (9s~k | f)PG t - k \F t - k (9t-k I f)df 



1 f°° 

-( 9s - k+9t - k )/2 / e -f Io (^g—f)I (^g—f)df. 



(40) 



1 f°° 



October 28, 2011 



DRAFT 



20 





(a) 



(b) 



(c) 



r -. > 




(d) 



(e) 



(f) 



Fig. 9. SSIM maps of (a) noisy version of slice #4 and (b)-(f) reconstruction results obtained using TVDN, WDN, GNLM, NLMS and NLMR, 
respectively. 



1 f°° 

_ e -(9s- k +a t - k )/2 ye -V> J Q {l y /g^y)j Q [ % Jg^y)dy, 



(41) 



where Jo is the zero-order Bessel function of the first kind, and i = y/^ T. Subsequently, using equation (2.32) of 
[421, we have 



SNL 



!) - L-(ff.-*+*-*)/4 Jn f yg«-fcgt-fc \ 



-(Ss-k+3t 



(42) 



where the last line stems from the fact that 7q(— x) = Iq(x). 



Appendix B 
Rational SM for Rician statistics 



In the case of Rician statistics, the rational similarity measure is given by 

SNL i 2 t,/c =PA a _ fc /A t _ fe |M s _ fe ,M t _ fc (l I m s _ fe ,m t _ fe ) = 

aPA s _ k \M 3 _ k (a | m s -k)pA t _ k \M t _ k {a I m t - k )da. 



(43) 
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(a) 



(b) 



(c) 




(d) (e) (f) 

Fig. 10. (a) Sagittal MRI scan; (b)-(f) Reconstruction results obtained using TVDN, WDN, GNLM. NLMS and NLMR, respectively. 



Proceeding with a similar derivation as in Appendix [A] we obtain 

/>oo 

SNL i 2 t,/c = / aPM a _ k \A^ k {m s -k | a)p Mt _ klAt _ k (m t -k I a)da 
Jo 



a" 



_ m s - k m t -k -( m ;_ t +mL t )/2^ / „„-a 2 /<r 2 r ( m s-ka \ T (m t - k a 



Substituting j/ = a/c in ( |44| results in 



SNL?i fc = m '-^-* e -( m ;->+^->)/^ /" ye -v* Io (™°=k y ) Io (™±± y ) dy . 



s,t,k 



(7 Z 



Finally, using formula (2.32) of (42], we have 

(2) 

S,t,k 2cr2 



II 



(44) 



(45) 



(46) 
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